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This paper presents techniques and results for simulations of unequal-mass, non-spinning binary 
black holes with pseudo-spectral methods. Specifically, we develop an efficient root-finding procedure 
to ensure the black hole initial data have the desired masses and spins, we extend the dual coordinate 
frame method and eccentricity removal to asymmetric binaries. Furthermore, we describe techniques 
to simulate mergers of unequal-mass black holes. The second part of the paper presents numerical 
simulations of non-spinning binary black holes with mass ratios 2, 3, 4 and 6, covering between 15 
and 22 orbits, merger and ringdown. We discuss the accuracy of these simulations, the evolution of 
the (initially zero) black hole spins, and the remnant black hole properties. 
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I. INTRODUCTION 

Numerical simulations of the inspiral and coalescence 
of two black holes [1] are an important tool for exploiting 
upcoming gravitational wave detectors such as Advanced 
LIGO, VIRGO, and LCGT/KAGRA gHB]. Increasingly 
larger sets of simulations have begun to explore the pa- 
rameter space of binary black holes (BBHs), most no- 
tably through the NINJA [7HS] and NRAR [TO] collabo- 
rations. 

One important subset of this parameter space com- 
prises non-spinning BBHs. Head-on collisions have been 
studied first [TTJ [12] , followed by simulations of inspiral 
and coalescence of binaries that start in a quasi-circular 
orbit. One well-studied phenomenon is the kick imparted 
to the remnant black hole as a result of the collision of 
unequal-mass black holes |T3l [14] ; the form of this kick as 
a function of the initial black hole masses is constrained 
by symmetry considerations [15 . Numerical simulations 
of non-spinning BBH systems also formed the basis of an- 
alytic waveform models and applications to gravitational 
wave data analysis [TSH2"2"] . tuning of effective-one-body 
waveform models [23H26] . multipolar analysis [57] EB] . 
and investigations into the periastron advance of binary 
black holes [291 123 ■ Recently, the range of mass ratios 
covered by unequal-mass binaries has been extended to 
mass ratios 10:1 [3JJ and up to 100:1 [52H54] , 

Numerical simulations are still too computationally 
expensive to include enough binary orbits for data 
analysis. Therefore, simulations are matched to post- 
Newtonian inspirals to obtain "hybrid" waveforms of suf- 
ficient length. This matching must be done early enough 
in the inspiral so that the post-Newtonian expressions 
are still accurate. During the last year, it has become 
increasingly apparent that current numerical simulations 
are still not long enough to provide an accurate match: 
the frequency range where post-Newtonian and numeri- 
cal waveforms are matched with each other is currently 
so high that neglected higher-order terms in even state- 
of-the-art post-Newtonian models lead to a noticeable 



impact on data analysis [TU I35H40] . 

Unfortunately, the computational expense of a BBH 
inspiral is a steep function of its initial frequency. For 
instance, at lowest post-Newtonian order [41], a BBH 
inspiral starting at an initial frequency Qi merges at a 
time 

T=^rT 1 (Mfi i )- 8/3 M ) (1) 

where M is the total mass of the binary and rj its sym- 
metric mass ratio r\ = M\Mij (M\ + M2) 2 . So even if the 
computational expense were proportional to the evolu- 
tion time T, it would be expensive to significantly reduce 
Qi] in practice the situation is even worse because the 
computational expense (for a given accuracy) increases 
superlinearly with T. Therefore, long numerical inspi- 
ral simulations (lasting > 10 orbits) are rare, and are 
generally available only for equal-mass binaries without 
spin [32], or with equal spin magnitudes parallel to the 
orbital angular momentum [43 [ 144] . 

This paper revisits simulations of non-spinning 
unequal-mass binary black holes, and describes accurate 
many-orbit waveforms, including subdominant (£, m) 
modes. Our simulations are performed with the Spectral 
Einstein Code SpEC [45 , a multi-domain pseudo-spectral 
evolution code. There are several motivations for this 
work. First, we present an efficient technique to perform 
10-dimensional root-finding that is necessary to construct 
BBH initial data with specified masses and spins. Sec- 
ond, we present algorithms for simulations of unequal- 
mass BBH systems with spectral methods. Third, we 
present and carefully discuss a series of long duration, 
high-accuracy, unequal-mass non-spinning BBH simula- 
tions, lasting between 15 and 22 orbits. These simu- 
lations extend the parameter space covered by spectral 
BBH evolutions, and improve in length and accuracy al- 
ready existing simulations which use alternative numeri- 
cal techniques. The simulations presented here also pro- 
vide additional data points for remnant masses, spins and 
kick velocities, which we compare with already published 
calculations and analytical models. Finally, we provide 
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a study of tidal spin-up of initially non-spinning black 
holes. 

This paper is organized as follows. Section |TT| presents 
details of our numerical implementation. First, the 
quasi-circular, quasi-equilibrium initial data |46l 147] re- 
quire root-finding to adjust free parameters so that af- 
ter the initial data construction, the black holes have 
specified masses and approximately zero spins - we in- 
troduce an efficient algorithm for performing this root- 
finding. Second, we extend the dual-frame approach [48] 
to unequal-mass binaries, and discuss how we choose or- 
bital parameters that result in inspirals of orbital ec- 
centricity e < 10~ 4 . Finally, we describe the handling 
of merger and ringdown, improving on previous treat- 
ments [12 HS1 HH] of black hole mergers performed with 
spectral multi-domain methods. Section [lll| presents nu- 
merical results for mass ratios 2, 3, 4, and 6. These in- 
clude results of convergence tests, discussion of the black 
hole spin, detailed analysis of the leading higher-order 
modes of the emitted gravitational waveform, and discus- 
sion of the properties of the remnant black hole: mass, 



spin and recoil velocity. Section IV summarizes and dis- 
cusses our main results. 

We note that the simulations presented here have al- 
ready been used in the following published work: fitting 
effective-one-body-models 25, 26, and measuring the pe- 
riastron advance for BBHs [53] . They have also been con- 
tributed to the Ninja2 [5] and NRAR projects [TU ] . Fur 
ther, the formalism for setting initial data (cf. Sec. 
and for eccentricity removal (cf. Sec. II D) was emp 
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II. FORMALISM & NUMERICAL METHODS 

A. Overview 

Our goal is to compute the last ~ 20 inspiral orbits, 
merger and ringdown of binary black holes with mass 
ratio q = Mi/M 2 > 1, negligible spins of the black holes, 
and vanishingly small orbital eccentricity. This requires 
a rather complex sequence of steps: 

1. Choose the physical black hole masses Mi, M 2 . 

2. Decide on the initial coordinate separation Dq, and 
choose tentative values for the orbital frequency 
S7o and its time derivative, parameterized by do = 
D(t)/D (for instance, based on post-Newtonian 
formulae) . 

3. Fine-tune the 10 parameters that enter the initial 
data so that the initial data contain black holes 
with desired masses, desired spins (here, zero), and 
vanishing center-of-mass motion. 

4. Perform a short evolution lasting 2-3 orbits of the 
resulting initial-data set. 



5. From the evolution in Step [4j extract information 
about the orbit of the binary and estimate the or- 
bital eccentricity e. If e is unacceptably large, cor- 
rect fio and do and go back to step [3] 

6. If the orbital eccentricity e is sufficiently small, con- 
tinue the evolution through the remaining inspiral 
(for the current paper, we require e < 10~ 4 ). 

7. Simulate plunge, merger and ringdown. 

In order to accomplish our goal, we needed to make sev- 
eral refinements to previous procedures used in SpEC for 
equal-mass [J5J 03 EH E2] and more generic (including 
q — 2 unequal-mass) [49] BBH simulations. These are: 
Step [3] was not necessary in previous evolutions of sim- 
pler configurations, and is explained in detail in Sec. II B 
below. Modifications to the inspiral evolutions in Step 4 
are detailed in Sec. |II C| Eccentricity removal in Step 5 
is generalized to mass ratios q ^ 1 in Sec. |IID| Improve- 
ments to the merger and ringdown phases (Step [7]) are 
described in Sec. II E| Finally, Sec. II F summarizes code 
infrastructure that has not changed since earlier simula- 
tions; examples are apparent horizon finders and wave 
extraction. 



B. Initial data 

Quasi-equilibrium binary black hole initial data [If)] 
|4"TI 155] are constructed with the conformal thin sandwich 
method E3] [55] . This formalism results in a set of five 
coupled non-linear elliptic equations, which are solved 
numerically with a multi-domain pseudo-spectral collo- 
cation method [SS] . 

As in earlier work, we employ the simplifying assump- 
tions of conformal flatness and maximal slicing. Thirteen 
further real parameters uniquely determine the complete 
initial data set. The orbital characteristics are deter- 
mined by the three parameters Dq (coordinate separa- 
tion), J7o (orbital frequency), and do (radial expansion 
factor) ; their choice will be discussed in detail in Sec 
The remaining 10 parameters 
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(r 1) r 2) fi 1 ,n 2 ,x,y) 



(2) 



are the radii ri,r 2 of the excision spheres, the angular 
velocities of the horizons, Qi, f2 2 , and the coordinate 
centers of the excision spheres, parameterized by X and 
Y via ci = (X, Y, 0) and c 2 = (X - D , Y, 0). We assume 
that the black holes start in the xy plane, with orbital 
angular frequency parallel to the z-axis, i.e. the vectorial 
orbital frequency is written as £7 = (0,0, Qq). 

The physical parameters (masses, spins, linear momen- 
tum) can only be computed after the constraint equa- 
tions are solved, whereas the initial data parameters u 
must be chosen beforehand. Therefore, 10-dimensional 
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root-finding is required, to satisfy 



F(u) ={Mx-M[, M 2 
=0. 



M' 2 , 

ADM) r ADMi 



(3) 



Here, Mi j2 , Xi, 2 , and Padm are, respectively, the masses, 
dimensionless spins, and total linear momentum, de- 
termined from the solution of the constraint equations, 
whereas M[ 2 and x'i 2 are the desired masses and dimen- 
sionless spins of the black holes. We also demand that 
the initial ADM linear momentum Padm vanish. The x— 
and y— components of Padm are controlled by the choice 
of Y and X, respectively. Its z -component Padm van- 
ishes by symmetry z — > — z (in generic spinning cases, 
this will no longer be the case). 

In this paper, we will evolve only non-spinning black 
holes such that x'i 2 = 0' but we present the root-finding 
for generic spins. 

Each function evaluation F_(u) requires solving the el- 
liptic constraint equations. At high resolutions, this re- 
quires a few hours of wall-clock time. Because root- 
finding with standard techniques such as the Newton- 
Raphson method |57j requires many function evaluations 
to compute the Jacobian, this would result in inconve- 
niently long run tirne^] To reduce computational ex- 
pense, we replace the exact Jacobian dF/du by an ap- 
proximation J a and perform a Newton-Raphson itera- 
tion employing J a- That is, given parameters u^ k \ im- 
proved parameters are determined by 



Au 



,(*) 



(4) 



where J a is evaluated at u( k \ 

Efficiency of this technique hinges crucially on the 
quality of the approximated Jacobian J a- We compute 
J A based on considerations that are valid for single black 
hole initial data, and/or Newtonian gravity. Specifically, 
for conformally flat single black hole initial data with 
maximal slicing, the mass is proportional to the radius 
of the excision sphere; therefore, we take 



QM A 
dr A 



Ma 

TA ' 



A = 1,2. 



(5a) 



Furthermore, for Kerr black holes with small spin, the 
dimensionless spin parameter \ is related to the angular 



1 In earlier work on equal-mass binaries with equal aligned spins, 
this root-finding was not performed. For those configurations, 
symmetry implies r\ = r<2, S^i = &2 and X = Y = 0. The 
radii ri = ri were chosen to be some fixed value, and the final 
black hole masses were simply measured (rather than controlled). 
For the non-spinning simulation |52| . f2i 2 were fixed at their 
values from quasi-circular non-spinning initial data |46| ; for the 
spinning simulation 1431 . f2i = Q2 was chosen parallel to the at- 
axia, and the resulting black hole spin was just measured (rather 
than controlled). 



frequency of the horizon fin by \ — 4Mfijj , where M is 
the mass of the Kerr black hole. For BBHs, the horizon 
frequency £Ia measures spin in addition to co-rotation, 
so that xa = 4Ma(S1a — Hq), from which follows 



dXA 
dr A 



XA 
T A ' 



dx/ 



A = 1,2. 



(5b) 



Finally, in Newtonian gravity, the linear momentum is 
given by P = MiQq X ci + M 2 f2o x ^2- Substituting in 
^0 = (0,0, Qo), ex = (X,Y,0), c 2 = (X-D ,Y,0), one 
finds 



Qpx 

dr\ 

Qpx 

liY 
gpv 

dr\ 
gpy 

~dX 



Mi 



n Y, 



3P X 



-{Mi 
M x 

(M 1 + M 2 )n 



QnX 



dr 2 

M 2 )n , 

gpy 



M2 
r 2 



dr 2 



M2 

T2 



-UnY, (5c) 
(5d) 

Q (X-D ), (5e) 
(5f) 



Equations (5a|-fl5f|) are the only non-zero components of 



J a- Because the Jacobian is so sparse, it is trivial to 
solve Eq. Q, and one obtains: 



Ar A = -r A 



M A - M' A 
M A ' 



An, 



AX 



Xa-x'a M a -M' a 

AM a 4M| 
_py 

r ADM 



A =1,2 (6a) 
Xa, A = 1,2 (6b) 



(M! + M 2 )n 

X(Mi-M{) + (X-D )(M 2 -M£) 



Mi +M 2 



— ADM 



(Mi + M 2 )n 

F(Mi - M[ + M 2 - M 2 ) 
+ Mi + M 2 



(6c) 



(6d) 



In these equations, primed quantities are the desired val- 
ues, whereas un-primed quantities are determined from 
the initial data computed from parameters yP^ . 

Fig. [T] demonstrates the efficiency of this procedure for 
two configurations. During the first iterations of root- 
finding, we solve the constraint equations only to lowest 
resolution. We begin to increase the resolution /ceii of the 
elliptic solver when the residual |P| falls within a factor 
of 10 4 of our target tolerance 10 -7 . Because solving the 
constraint equations at low resolution is very quick, the 
overall cost of the root-finding is dominated entirely by 
the solutions of the constraint equations at highest reso- 
lution, and thus, the entire root-finding adds only a small 
amount of wall-clock time. 

As is apparent in Fig. [l] the quadratic convergence 
of Newton-Raphson algorithm is lost because of the ap- 
proximations entering J a- We find roughly linear conver- 
gence where each iteration reduces the error by a certain 
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lation of Einstein's equations is fixed via a freely speci- 
fiable gauge source function H a that satisfies the con- 
straint 

= C a = T ab b + H a , (7) 

where T a f, c are the spacetime Christoffel symbols. During 
the inspiral, we choose H a as in Refs. [4"2l 14"51 152] . 

In order to treat moving holes using a fixed grid, we 
employ multiple coordinate frames [48]: the equations 
are solved in an 'inertial frame' that is asymptotically 
Minkowski, but the grid is fixed in a 'grid frame' in which 
the black holes do not move. The motion of the holes 
is accounted for by dynamically adjusting the coordinate 
mapping between the two frame^] This coordinate map- 
ping differs from our earlier work, and is described below 
in Sec. [HUH 

Furthermore, the choice of constraint damping pa- 
rameters is important for stability, and it discussed in 

Sec. Etna 



FIG. 1. Residual \F\ of the initial-data root-finding proce- 
dure (top panel) vs. iteration number. The lower panel indi- 
cates the numerical resolution of each elliptic constraint solve, 
with 5 being highest resolution (each increase of this integer 
corresponds to adding a certain number of basis-functions, cf. 
Fig. 6 of 07]). 



factor. The convergence rate depends on how closely J a 
resembles the exact Jacobian. Convergence is not ex- 
actly linear, because we delay increasing the resolution 
of the elliptic solver until as high k as possible, to gain 
maximum speed-up from the lower resolution solutions. 



C. Evolution of inspiral phase 

The Einstein evolution equations are solved with 
the pseudo-spectral evolution code (SpEC) described in 
Ref. [H|. This code evolves a first-order representa- 
tion [58] of the generalized harmonic system |59H61j and 
includes terms that damp away small constraint viola- 
tions 58, EU [62]. The computational domain extends 
from excision boundaries located just inside each appar- 
ent horizon to some large radius, and is divided into sub- 
domains with simple shapes (e.g. spherical shells, cubes, 
cylinders). No boundary conditions are needed or im- 
posed at the excision boundaries, because all character- 
istic fields of the system are outgoing (into the black hole) 
there. The boundary conditions on the outer bound- 
ary [SSI IM1 El] are designed to prevent the influx of 
unphysical constraint violations [65 71 and undesired in- 
coming gravitational radiation [72, 73 , while allowing the 
outgoing gravitational radiation to pass freely through 
the boundary. Interdomain boundary conditions are en- 
forced with a penalty method [7U H5] ■ 

The gauge freedom in the generalized harmonic formu- 



1. Dual-frames & Control system 

SpEC utilizes two coordinate systems [48]: grid coor- 
dinates x l , in which the domain decomposition is fixed, 
and inertial coordinates x z , in which the black holes orbit 
around each other. The mapping between these coordi- 
nate systems is chosen such that in grid coordinates, the 
black holes remain centered on the excision spheres. In 
earlier simulations of equal-mass binaries [48, 5TJ[52], this 
map was chosen to be a rotation and an overall scaling. 
Unequal-mass binaries will acquire a kick in the orbital 
plane; therefore, we add a translation to the mapping 
between inertial and grid coordinates: 



x l = a(t)R\x l + T 



(8) 



Here, a(t) is the overall scale factor, T l = (T 2 ,Tf,0) 
represents the translation, and 



R 1 



o 1 



R,. 



cos cp — sin < 
sin <b cos ( 



(9) 



is the rotation matrix for a rotation by the angle <p(t) 
about the z-axis. The rotation and translation act only 
on the x— and y— coordinates, because a non-spinning 
unequal-mass binary is, by symmetry, confined to remain 
in the xy-plantj^ 

The mapping Eq. ([8]) is determined by four free- 
functions, X a ee {a(i),0(t),r 2 (t),Tf(i)} (where a labels 



2 All coordinate quantities (e.g. trajectories, waveform extraction 
radii) in this paper are given with respect to the inertial frame 
unless noted otherwise. 

s Spinning, unequal-mass binaries with both black hole spins par- 
allel to the orbital angular momentum will also remain in a fixed 
orbital plane. Our discussion applies equally well to these sys- 
tems. 
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the four functions). The functions X a (t) must be cho- 
sen dynamically such that the black hole horizons re- 
main centered on the excision boundaries. As described 
in Ref . [25] , this is accomplished through a control system 
that constantly monitors the location of the black holes, 
and dynamically changes the functions X a (t) appropri- 
ately. Such a control system is formulated most easily 
in terms of control parameters Q a = {Q a ,Q<j),Qx,Qy} 
which have the properties (i) that Q a = when the black 
holes are at their desired locations, and (ii) for small val- 
ues of Q a , changing the mapping-parameters X a changes 
the control parameters Q a according to 



dQa 
dX B 



-S a 



|8> 



for \Q a \ < 1. 



(10) 



The control parameters must be given in terms of the 
moving coordinates of the centers of the apparent hori- 
zons, c\ 2 , and they must vanish when c\ 2 are at the 
desired locations, namely, when they are at their va lues 
in the initial data (c\ 2) t _ Q - The derivatives in Eq. (10) 
are to be taken at constant inertial coordinates of the 
centers of the horizons. 
To begin, we define 



(A x (i), A y (t), A z (t)) = c^t) - c 2 (t), 

D(i)=[A»(t) + A*(i)] 



1/2 



(11) 

(12) 



Because of symmetries, A 2 is always zero, and will not 
be used. The control parameters for the expansion factor 
a(t) and the rotation angle <fi are given by 



V0 ~ D(t) ' 



(13a) 
(13b) 



It is straightforward to verify that Q a and Q</> satisfy 



Eq. (10) 



The control parameters for the translation are some- 
what more involved. We use the ansatz 



Q ;?' 
Qy 



= a(i)R, 



!JB 



M 



(13c) 



where M is a constant 2x2 matrix, and we demand that 



M commutes with R 



'0(t)' 



Because M and R^,(t) com- 



mute, Eq. ( 13c ) can be rewritten in inertial coordinates 
as 



Qx 

Qy 




(14) 



which makes it obvious that Q x and Q y satisfy Eq. (10). 

To close this discussion, we must compute the matrix 
M. The requirements that M commute with R^ and 
that Q x = Q y — for c\ 2 = (cl 2 ) t _ determine M 
uniquely: 



M = 



1 

Do 



XA,0 
VA,0 



~VAfi 
XA,0 



(15) 



The mapping given in Eq. (|Sj) and the control param- 
eters, given in Eqs. (13), are then combined with the 
feedback control system described in Ref. [JSJ m order to 
evolve the unequal-mass BBH through the inspiral phase. 



2. Constraint Damping 

In order to suppress violations of the generalized har- 
monic gauge constraint Eq. ([7j (cf. Refs. [SH[7E|), an d 
of the auxiliary constraints that arise from the reduction 
of the generalized harmonic evolution system to first or- 
der form (cf . Ref. [53 [77] ) , we introduce so-called con- 
straint damping terms in the generalized harmonic evo- 
lution equations (see [55]). These terms are proportional 
to the constraint damping parameters 70 and 72 . 

Simulations with mass ratios q = {2,3} were found 
to be stable with the same constraint damping parame- 
ters as those used in Ref. [S2]. However, for the higher 
mass ratios q = {4,6}, we encountered constraint viola- 
tions that grew exponentially on time scales of several 
100M. We found that toward the outer edges of the 
cylindrical subdomains, the constraint damping param- 
eters must be sufficiently large in order to suppress ex- 
ponential constraint growth. In the overlap between the 
inner spherical shells and the cylinders, an instability de- 
velops unless the constraint damping is sufficiently small. 
Furthermore, we were not able to achieve stable evolu- 
tions with 70 =72- After considerable experimentation, 
we settled on a sum of Gaussians: 

M 7o =8e-^/ 1 ' 3M ) 2 + 16e~^/ M > 2 + / fai _ field (r) (16) 
M 72 =8e"^/ 1 - 3M ) 2 + 40e-^/ M > 2 + / fax _ fie i d (r) (17) 

with far-field terms /far-field = 0.2 e -(''/ 60M ) 2 + 0.001. 
Here r\ and r 2 are the coordinate distances from the 
centers of each hole, and r is the distance from the origin. 
The choices Eqs. ( 16 1 and ( fl7| ) were found to work well 
even for q — {2, 3}, and all simulations presented here 
use them. 

We infer from these results that the domain decom- 
position with spheres overlapping cylinders is not always 
stable, and that stability depends sensitively on certain 
geometric details. Recent shorter simulations that do not 
have overlapping subdomains do not show such sensitiv- 
ity. However, the domain decomposition of spheres and 
cylinders is computationally more efficient, and therefore 
we employ it during long inspiral simulations. 



D. Eccentricity removal 

The procedure for eccentricity removal developed in 
Refs. [51] [52] assumed an equal-mass binary. Generaliza- 
tion to unequal-mass binaries is straightforward. As in 
Ref. [S5], we fit the radial velocity (represented by the 
time derivative of the proper separation s(t) between the 
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FIG. 2. Eccentricity removal for mass ratio q = 4. The 
main panel shows proper separation as a function of time, the 
inset radial velocity ds/dt. Initial data parameters based on 
TaylorT3 post-Newtonian approximation result in the black 
dotted line with e ~ 0.008. The red-dashed and the solid- 
blue lines represent two iterations of eccentricity removal, for 
a final eccentricity of e w 4 x 10 -5 . 



The orbital eccentricity is given by 

B 

SqUJ 7 



S-ds/dt 



(22) 



which is the same formula as for the equal-mass case. 

Overall, eccentricity removal works as well here as for 
the equal-mass cases considered previously. Fig. [2] shows 
that with each iteration, e drops by about a factor of 10. 
The most important factor for effective eccentricity re- 
moval is the quality of the fit. The fitting interval [ii,i2] 
can start only after transients due to junk radiation have 
decayed. However, because the fit is used to infer radial 
velocity and acceleration at time t = 0, the fitting in- 
terval needs to be sufficiently early in the run to allow 
accurate extrapolation from the fitting interval back to 
t = 0. Finally, the fitting interval needs to be long enough 
to allow a reliable fit of the frequency u), i.e. it needs to 
be longer than one period of the radial oscillations. In- 
clusion of the term quadratic in t in Eq. ( 19 ) significantly 



improves the quality of the fits and the effectiveness of 
the eccentricity removal. For the runs described here, we 
choose t\ on the order of 100M and £2 on the order of 
1000M. 



E. Evolution of merger & ringdown 



horizons) by the functional form 
ds 

-j-=v insp (t) + Bcos(ut + (t)). (18) 

Here Vi nsp (t) is a monotonic function varying on the 
(long) inspiral time scale; this function captures the 
desired zero-eccentricity inspiral driven by radiation- 
reaction. We take here the functional form 

Vinsp(£) = V + V\t + V 2 t 2 , (19) 

with three fitting parameters vq,Vi, v%. However, in more 
recent work [50] , we describe fitting functions that result 
in more robust behavior. The oscillating piece B cos(a;£-l- 
</>) captures superposed oscillation due to non-zero orbital 
eccentricity - the goal is to reduce the amplitude of this 
piece. 

For unequal masses, the black holes have different sep- 
arations from the origin, and therefore have different ra- 
dial velocities. To avoid dealing with each black hole 
independently, we consider the initial data specified in 
terms of a Hubble-like radial expansion factor do, which 
induces radial velocities proportional to the distance to 
the origin, v l r = clqx 1 at a coordinate location x l . The 
updating formulas become 

fionew = ^0 + 7T~ Sm (<? !) ): ( 20 ) 

2s 

donow = do cos(<£). (21) 

so 



The evolution algorithm for the inspiral described in 
Section |II C| fails when the black holes approach each 
other too closely This failure is caused by several factors. 
First, the gauge fields H a are chosen during inspiral to 
be time-independent in the grid frame. This works well 
for the inspiral because the solution (in the grid frame) 
is roughly time-independent near the black holes. Near 
merger, however, this gauge leads to the formation of 
coordinate singularities. Second, during inspiral, the ex- 
cision boundaries of the grid remain spherical, and do not 
change shape even though the individual apparent hori- 
zons become distorted as the holes approach each other. 
As the distortion of the apparent horizons increases, the 
mismatch between the excision boundaries and the ap- 
parent horizons eventually leads to a violation of the ex- 
cision condition, i.e., the condition that all characteristic 
fields of the hyperbolic system are outgoing (i.e. into the 
hole) at each excision boundary. Third, the overlapping 
domain decomposition used during the inspiral is prone 
to weak instabilities that cause no trouble during the in- 
spiral but drive rapidly growing modes after the solution 
becomes highly dynamical. 

To address these problems, we stop the simulation 
about 1.5 orbits before merger, and restart with a mod- 
ified algorithm. We change smoothly to a damped har- 
monic gauge [IHl [751 that slows down the formation 
of coordinate singularities. We also dynamically modify 
the coordinate mapping between the grid frame and the 
inertial frame so that the excision boundaries conform 
to the shapes of the apparent horizons [121 [35]. Fur- 
thermore, by monitoring the characteristic speeds of the 
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FIG. 3. Domain decomposition used for the plunge and 
merger for mass ratio q = 2. The thick blue lines represent 
subdomain boundaries in the z~0 plane. The region z > is 
not shown. Also not shown is the additional deformation of 
the grid near the black holes that matches the shape of the 
excision spheres to the apparent horizons. 



system, we dynamically vary the velocity (with respect 
to the horizon) of each excision boundary so as to en- 
sure that the characteristic fields are outgoing at these 
boundaries for all times; this characteristic speed control 
is also crucial for evolving BBHs with large spins [44] . Fi- 
nally, we run the simulation on a set of non-overlapping 
sub-domains consisting of topological cubes, cylindrical 
shells, and spherical shells. This domain decomposition 
is shown in Fig. [3] Each subdomain is distorted by a 
coordinate mapping so that the subdomains do not over- 
lap and so that the union of these subdomains covers the 
entire 3-dimensional region (minus two excised holes) in- 
side a spherical outer boundary RBdry of ord er a few 
hundred M from the source (see Section III C 2 where we 



compare runs with different values of RBdry)- More de- 
tails about the merger domain decomposition are given 
in the Appendix. It avoids certain instabilities that ap- 
pear for domain decompositions with overlapping grid 
close to merger jj^]. In addition, we choose a slightly 
higher resolution for the non-overlapping grid than for 
the overlapping grid used during inspiral, because the 
merger has features with a shorter length scale than in 
the inspiral. After the binary has reached about t ~ 2M 
before merger, we increase the resolution one last time, 
particularly in the region between the two holes[^] 

After a common apparent horizon forms, we regrid 
onto a new set of subdomains consisting of nested dis- 
torted spherical shells. The innermost boundary is just 
inside the common apparent horizon, and conforms to 
its shape. The outermost boundary is the same RBdry 



used in the merger. The matching of the ringdown to 
the inspiral is discussed in |49j . 



F. Relation to other SpEC simulations 

Several other SpEC simulations of binary black holes 
have been presented in the literature [101 |4"2"H141 152"] . In 
this section we briefly describe some computational de- 
tails common to all SpEC simulations, and we describe 
how some of the new computational infrastructure pre- 
sented here relates to these other simulations. 

Our apparent horizon finder expands the radius of the 
apparent horizon as a series in spherical harmonics up 
to some order L. We utilize the fast flow methods de- 
veloped by Gundlach JSU] to determine the expansion 
coefficients. The quasi-local spin S of each black hole 
is computed with the spin diagnostics described in [81] . 
We compute the spin from an angular momentum sur- 
face integral [S3] using approximate Killing vectors 
of the apparent horizons, as described in [81] [84] (see 
also [551 [BE])- We define the dimensionless spin by 



X 



_S_ 
M 2 ' 



(23) 



4 The processes of regridding, changing resolution, and changing 
the coordinate mapping have since been automated; this will be 
described in a future work. 



We extract gravitational waves from our simulations 
by two independent methods. We compute the Newman- 
Penrose scalar '5 4 using the same procedure as described 
in [5"H [52] . This involves constructing the correct con- 
traction of the Weyl curvature tensor at several finite- 
radius coordinate-spheres far from the source and pro- 
jecting into spin- weighted spherical harmonics. We also 
extract the Regge-Wheeler-Zerilli (RWZ) [B2 [88] gravi- 
tational wave strain hg m as formulated in Ref. [85]. The 
implementation of this formulation in the SpEC code is de- 
scribed in [90] (see also [26] and the appendix of [25] for 
further details). Both the ^4 and the RWZ waveforms, 
which are extracted at a series of finite-radius coordinate 
spheres, are extrapolated to infinite distance from the 
source [91] . The 'J 4 waveforms generally agree well with 
the (second time derivative of the) RWZ h^ m waveforms, 
although for some purposes RWZ is a better choice than 
^4 or vice versa. For example, computing strain from ^4 
requires two time integrations and careful choice of inte- 
gration constants, so it is simpler and less error-prone to 
instead use RWZ to compute strain. Similarly, comput- 
ing the recoil velocity requires either a time derivative of 
h£ m or a time integral of ^4; the time derivative amplifies 
noise in the waveform, and this affects the recoil velocity 
enough that it is better to use a time integral of $4 for 
that purpose. 

In parallel to the present work, superposed Kerr-Schild 
initial data[HU [H21 [53] have been developed and applied 
to SpEC simulations of black holes with high spins [4T)ll44| . 
The algorithmic improvements discussed in the present 
work are generally compatible with superposed Kerr- 
Schild simulations. Specifically, the root-finding proce- 



dure discussed in Sec. II B can be applied to superposed 



Kerr-Schild initial data. This requires a change of free 
parameters from excision sphere radii to masses of the 
conformal black holes in the superposed Kerr-Schild ini- 
tial data. Early tests indicate that the root-finding pro- 
cedure works satisfactorily. However, more exhaustive 
tests, especially for high spin systems, will be necessary. 

The control system discussed in Sec. |II C 1| is applica- 
ble to any non-precessing simulation, independent of the 
type of initial data. The choice of gauge source functions 
H a (equal to the values in the initial data, with appropri- 
ate coordinate transformations applied [121 H31 [55] ) does 
not work for simulations with moderate or large spins; 
such simulations use active gauge conditions already dur- 
ing the inspiral, see e.g. [3Dj. Furthermore, moderate to 
high spin simulations require use of a non-overlapping do- 
main decomposition during the inspiral to avoid certain 
grid instabilities. We have no reason to believe that the 
more complex and computationally more expensive tech- 
nology for high-spin systems might fail for the present 
non-spinning simulations. We have not tested this, be- 
cause the methods presented here are more efficient for 
the systems being studied here. Techniques for handling 
the merger, as described in Sec. |IIE| and the Appendix, 
are common between the high-spin simulations and the 
simulations presented here. 
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III. RESULTS 



Overview 



In this section, we present the results of our simula- 
tions of non-spinning binary black holes with mass ratios 
q = 2,3,4,6. These simulations contain long inspirals 
(15 to 22 orbits), merger, and ringdown. To achieve 
our desired number of inspiral orbits, we compute the 
initial coordinate separation D using Taylor T3 post- 
Newtonian predictions [JT], an d then proceed to the ec- 
centricity removal procedure as explained in Sec. II D 



Our final parameters for the initial data set are summa- 
rized in Table [TJ and Fig. [4] shows the trajectories of all 
our runs through inspiral, the formation of a common 
apparent horizon, and merger. 



FIG. 4. Orbital trajectories for mass ratio q = 2, 3, 4, 6. For 
all mass ratios, the trajectory of the larger hole is represented 
by a dashed blue line, and that of the smaller hole by a solid 
red line. 



is the relative change of M- m (t). Convergence is clearly 
apparent, and the irreducible mass is constant to within 
a few parts in 10 6 at the highest resolution, except im- 
mediately before merger. During the first ~ 100M, the 
black hole mass increases by about 1 x 10~ 6 . Since this is 
below the numerical error during inspiral shown in Fig. [5] 
we define our mass scale by 



M = M irr (0) 



(25) 



for all mass ratios. 



B. Mass calibration 

A mass scale M by which all data are rescaled is de- 
fined as follows. Consider the sum of the two irreducible 
masses, defined from the areas Aah i and A^h i of the 
apparent horizons, 

Root-finding during construction of the initial data en- 
sures Af irr (0) = 1. Figure [5] presents convergence data 
for the irreducible mass during the simulations. Plotted 



C. Accuracy 

1. Phase convergence 

One of the goals of the present work is to calculate long, 
accurate waveforms for the dominant and top subdomi- 
nant gravitational wave modes - (2,2), (3,3), and (2, 1) 
- from unequal-mass binary black hole simulations. The 
top subdominant modes are those with the largest peak 
strain amplitude. To determine the accuracy of these 
waveforms, we perform convergence studies of RWZ-ft^ m 
at a particular extraction radius. 
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3 


31 
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20.3077 


-66.08 
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402M 


400M 


4 


31 


0.97792(2) 


0.47160(10) 


157(2) 


6 


19.35244 
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13.0000 


0.995968 


0.5157 


572M 


569M 


4 


43 


0.98547(5) 


0.37245(10) 
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TABLE I. Runs considered in this paper, with q — 1 from Ref. [42] included for completeness. Initial data parameters are 
orbital frequency fio, the expansion factor do, and the coordinate distance between the black hole centers Do- Furthermore, the 
initial and final radii of the outer boundary are given (7?Bdry is decreasing during the evolution, cf. [52]), as well as the initial 
orbital eccentricity Eds/dt and the number of gravitational wave (GW) cycles before the peak of I/122I, Ngw- The last three 
columns denote the Christodoulou mass, dimensionless spin, and kick velocity of the merged black hole at the end of ringdown. 




FIG. 5. Convergence of the irreducible mass. Plotted is FIG. 6. Phase convergence of RWZ-h (2,2) modes for inspiral- 
5M irr = [M iTT (t) ~ M irr (0)] /M irr (0) for three different numer- merger-ringdown waveforms. Shown are the phase differences 
ical resolutions (TV = 3 4 5) between a given resolution and the highest resolution. 



All simulations are run at three different resolutions, 
labeled N — 3,4,5. For all three resolutions, the 
RWZ gravitational waveforms at a finite extraction ra- 
dius (.fl^rt = 338M for q = 2,3,4 and H oxt = 460M 
for q = 6) are computed. We decompose the complex 
spherical harmonic modes into real- valued amplitude and 
phase: 

h lm {t) = A lm (t)exp(i(j) lm (t)). (26) 

We next compute differences A<f>i m (t) between different 
resolutions without any time shifts, 

A^'(t) = C(*)-^(*)= (27) 

where the superscripts N and N' refer to the numerical 
resolutions being considered. Finally, for ease of presen- 



tation, we time-shift the phase differences to align con- 
vergence tests of different mass ratios at their respective 
times of peak amplitude of the /122 mode, i pea k22- 

Phase differences for the dominant (2, 2) mode are plot- 
ted in Fig. [6] Note that this figure shows only the part of 
the simulation around merger time. During the earlier in- 
spiral, the phase errors are lower. It is apparent from this 
plot that the phase accuracy deteriorates with increased 
mass ratio, albeit quite slowly. This is expected, as simu- 
lations become numerically more difficult with increased 
mass ratio, owing to the smaller gravitational wave (GW) 
flux, and the smaller length scale of the small black hole. 
Nevertheless, the phase accuracies of all the new simula- 
tions presented in this paper are comparable to that of 
the equal-mass, zero spin simulation presented in Scheel 
et al [32] j with the simulations at low mass ratios (q = 2) 
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being somewhat more accurate, and those at higher mass 
ratios (q — 3, 4, 6) somewhat less accurate. 

Note that during merger and ringdown, the three res- 
olutions of the q = 2 simulation do not follow the usual 
pattern indicating convergence. There are a few possible 
reasons for this. One is that for q = 2, the truncation er- 
ror as a function of resolution may change sign near one 
of the resolutions N = 3, 4, or 5, thus producing an arti- 
ficially small truncation error and skewing the test shown 
in Fig. [6] Another possibility is that the unusual pattern 
is caused by small differences in gauge or domain decom- 
position between different resolutions: as explained in 
Section pi E| we change the gauge and domain decompo- 
sition about 1.5 orbits before merger, but these changes 
occur at slightly different times for different resolutions, 
and this time offset will introduce a small non-convergent 
error. Note also that the q = 2 case appears to have a fac- 
tor of three smaller truncation error than any previous 
long SpEC simulation, so this case may reveal small error 
sources that may not have been evident in previous simu- 
lations. Figure [6] shows a feature in the q — 6 simulation 
around t ~ — 180M. This arises because the phase dif- 
ference between N = 4 and N = 5 simulations changes 
sign. 

Convergence tests for the two leading subdominant 
modes (2, 1) and (3, 3) are presented in Fig. [7] Dur- 
ing the inspiral, the phase errors of the (2,1) mode are 
approximately half as large as those for the (2,2) mode, 
whereas the errors in the (3,3) mode are approximately a 
factor 1.5 larger. This scaling is reasonable, as all three 
GW modes are determined primarily by the orbital phase 
evolution. The gravitational wave mode (l,m) proceeds 
through m cycles for each orbit; hence, the GW phase 
errors of different modes should be proportional to m. 
During merger and ringdown, the observed phase errors 
behave differently: A033 is larger than A022 for all mass 
ratios, whereas A02i is similar in amplitude to A</>22- 
Fig. [7] shows noise in the (2,1) convergence test, starting 
about 150M before peak amplitude. Presumably, the 
noise in the phase is more prominent in the (2,1) mode 
because of the small amplitude of this mode. 



2. Effect of location of outer boundary 

The simulations presented here are of such long dura- 
tion that the black holes are in causal contact with the 
outer boundary for a large portion of the evolution. The 
question therefore arises: are the results affected by our 
choice of outer boundary conditions? Ideally, the grav- 
itational waveforms computed on a truncated computa- 
tional domain with an artificial outer boundary should 
not have errors introduced by the boundary conditions 
themselves - either from spurious reflections of gravita- 
tional radiation or from constraint violations at the outer 
boundary. The extent to which this is achieved indicates 
the degree to which the outer boundaries are "absorb- 
ing" (see e.g. Refs. [Ml HI ES |SD] ) . The outer bound- 




FIG. 7. As Fig. [6] but for the subdominant modes (3,3) and 
(2,1). 



ary conditions used in our simulations are (i) constraint- 
preserving and (ii) freeze the Weyl scalar ^0 to its ini- 
tial value. These "semi-absorbing" boundary conditions 
are the simplest in a hierarchy of increasingly absorb- 
ing boundary conditions, described in detail in Sec. 4.2 
of [72]. 

To evaluate the impact of the artificial outer boundary 
on our simulations, we repeat the N — 4 simulations for 
each mass ratio with two additional outer boundary radii, 
^ciosc an d -Rfar, where the distance to the outer bound- 
ary is changed only by adding or removing outer spherical 
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FIG. 8. Effect of the outer boundary location. Shown 
are phase differences of h'22 between simulations with outer 
boundary radii given in Table [TTJ The solid lines give the dif- 
ferences between "normal" and "far" boundaries; the dashed 
lines give the differences between "close" and "far" bound- 
aries. For clarity, q = 2, 3, 4 are offset vertically by multiples 
of 0.02rad and q — 2, 3 are offset horizontally by multiples of 
300A/. 7? denotes the extraction radius of each simulation. 



shells in our domain decomposition. The different outer 
boundary radii are listed in Table [TTJ The /122 waveforms 
are extracted from these simulations, and phase differ- 
ences between runs with different outer boundary radii 
are computed and plotted in Fig. [8j The plotted phase 
differences are oscillatory during inspiral, indicating that 
the runs being compared have slightly different orbital 
eccentricities. Around merger, a systematic phase differ- 
ence appears of a few times O.Olrad for the near bound- 
ary and < 0.005rad for the normal boundary location. 
During ringdown, the gravitational wave amplitude de- 
cays exponentially and the calculation of the phase be- 
comes increasingly noisy. We truncate the plotted data 
when the amplitudes of the waves have decayed to 1% 
of their peak values. It is evident from Fig. [8] that the 
transparency of the outer boundary diminishes as the 
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TABLE II. Radii of the outer boundary for the runs with 
different outer boundary locations (in units of the initial sep- 
aration Do)- Also given is the ratio Q of spurious reflections 
from the (2,1) mode relative to those from the (2,2) mode, cf. 
Eq. (29). 



distance to the boundary decreases. For our "normal" 
boundary radius, the phase error due to the boundary is 
< 0.005rad (when compared to the far location), which 
is negligible relative to the truncation error presented in 
Fig. [6] On the other hand, moving the boundary from the 
"normal" to the "close" location increases phase errors 5 
to 10 times. 

We can relate the phase errors reported in Fig. [8] to 
the expected reflection coefficients of our semi-absorbing 
boundary conditions as analyzed in Ref. [72] ■ The 
quadrupolar wave (£ — 2) reflection coefficient 02 for 
freczing-4'o plus constraint preserving boundary condi- 
tions is given by Eq. (89) of Ref. [72] • In the limit of large 
boundary radius /ci?Bdry 1 (where k is the wavenumber 
of the outgoing wave), the reflection coefficient reduces 
to 



a 2 = ^ (fc-RBdry) 4 • 



(28) 



"Near" boundaries are a factor ~ 1.6 closer than "nor- 
mal" boundaries; therefore, the reflection coefficient will 
be larger by a factor 1.6 4 ~ 6.5, consistent with the ob- 
served increase of phase errors by a factor 5-10 in FigJH] 
Moreover, according to an argument given in Ref. [94J . 
the phase error due to reflection of the (2, 2) mode of 
the outgoing radiation should be roughly equal to a 2 
times the total accumulated phas^J For the q = 2,3,4 
simulations with normal boundary locations, we have 
fci?Bdry ~ 18 and cr 2 ~ 1-3 x 10~ 5 . The - 30 GW- 
cycles of inspiral correspond to 2 2 ~ 200rad, so that 
02^22 ~ 0.003rad, in broad agreement with Fig. [8] 

For unequal-mass BBHs, it is important to consider 
reflection coefficients for higher-order modes, since the 
amplitude of these modes relative to the dominant (2,2) 
mode increases with mass ratio (see Fig. 10). For ex- 



ample, the reflection coefficients for both the (2,1) mode 
and the (2,2) mode are given by Eq. ( [28) , but the (2,1) 
mode has twice the wavelength of the (2,2) mode, reduc- 
ing fci?Bdry by a corresponding factor of 2. Consequently 
the reflection coefficient 021 of the (2,1) mode is a factor 
2 4 = 16 times larger than the reflection coefficient 022 
of the (2,2) mode. If we assume that the impact on the 
phase error is proportional to the amplitude of the re- 
flected waves, then the relative importance of reflections 
of the (2,1) mode and the (2,2) mode is given by the ratio 



Q 



m—l.m—2 



^21^21 
^22022 



(29) 



where A21 and A22 are the amplitudes of the (2,1) and 
(2,2) modes, respectively. Note that in the limit of large 
radii, Q m =i,m=2 is independent of boundary radius (be- 
cause i?Bdry cancels out of the ratio (721/1722) and inde- 
pendent of GW extraction radius (because the extraction 



5 Depending on assumptions, <T2 may be raised to a power close to 
unity, cf. Eq. (17) of Ref. [94]. 
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FIG. 9. Gravitational waveforms for q = 1, 2, 3, 4, 6 that have been extrapolated to infinity. The black curve is the /t22 mode, 
the red line is the /133 mode, and the blue line is the /121 mode. Only the real parts have been plotted. The a;-axis has been 
time-shifted so that indicates the merger, as determined by the peak of the extrapolated /122 mode waveform, for each mass 
ratio. The left-hand panels show the full coalescence: inspiral, merger, ringdown. The right-hand panels show a close-up of the 
merger and ringdown. Only the /t22 mode is shown for q = 1, since the odd-m modes do not appear here. 



radius cancels out of the ratio ^21/^22 )• Looking up the 
amplitudes of the (2,1) and (2,2) modes from Fig. [ToJ and 
using CT21 /022 = 16 results in the numerical values shown 
in Table [TT] (note that for these calculations, the ampli- 
tudes were taken at a specific time during the inspiral 
when they are still fairly constant). From this table, we 
conclude that with our semi-absorbing (constraint pre- 
serving plus freezing- \l/o) boundary conditions, the im- 
pact of the (2,1) reflections on the overall phase error 
is comparable to that of the (2,2) reflections, especially 
as the mass ratio increases to q — 4 or higher. With 
boundary conditions that are less than semi-absorbing, 
the error contributions would be even higher. 



D. Properties of gravitational radiation 

Fig. [9] shows the waveforms for our 15-orbit inspiral, 
merger and ringdown, as measured by (R/M)hi m . All 
these waves have been extrapolated to infinity. We show 
the top three modes: (2, 2), (3, 3), (2, 1). Notice that the 



amplitude of the (2, 2) mode decreases as the mass ra- 
tio increases, but the amplitudes of the other modes stay 
approximately the same. Further notice that the wave- 
length of the (2, 1) mode is about twice that of the (2, 2) 
mode. This is a general property: for a given £, the 
wavelength of the waveform is typically proportional to 
l/\m\. 



The relative importance of the (3, 3) and (2, 1) mode 
amplitudes to that of the (2, 2) mode is shown for the 
inspiral and merger in Fig. [To] (top panel: (3, 3) mode, 
bottom panel: (2,1) mode). This figure clearly shows 
that the higher order modes grow in relative significance 
as the mass ratio increases. At frequency Mio 22 = 0.06, 
the ratio A 33 /A 22 ranges from 0.08 (for q = 2) to 0.16 
(for q = 6), and A 21 /A 22 from 0.04 (for q = 2) to 0.08 
(for q = 6). At the peak of the h 22 waveform (indicated 
by the filled circles in Fig.llO]), A 33 /A 22 = 0.14 for q = 2 
and 0.28 for q = 6; ^21/^22 = 0.09 for q = 2 and 0.20 
for q — 6. 
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less gauge dependent) via M/b = (Afil) 2 / 3 , one obtains 
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FIG. 10. Amplitudes of /133 (top) and /121 (bottom) modes, 
normalized by the amplitude of the leading /i22 mode, for mass 
ratios q = 2, 3, 4, 6. Relative amplitudes are plotted versus 
the frequency of the /122 mode. The filled circles indicate the 
frequencies where the amplitude of the /122 mode peaks. 



E. Black hole Spin & Tidal spin-up 



We measure black hole spins by a surface integral on 
the apparent horizon that utilizes approximate Killing 
vectors computed from a minimization principle [81] . We 
denote the dimensionless spin by \A = Sa/M\ where 
A=l indicates the more massive black hole, and A = 2 
the less massive one. At t = 0, both black hole spins 
are very small: Xi(t = 0) < 1CP 8 . This is expected since 
Xa = is enforced as part of the initial data construc- 
tion, cf. Sec. |II B| During the initial relaxation of the 
initial data, the black hole spins increase to a few parts 
in 1CP 7 . Subsequently, \i slowly increases during the 
inspiral (with spin rotation axis parallel to the orbital 
angular momentum). This increase is convergently re- 
solved, as shown in the left panel of Fig.[TT] In contrast, 
the spin of the smaller black hole xi remains closer to 
zero, as shown in the right panel if Fig. 11 For mass 
ratios q = 3,4,6, xi is consistent with zero within trun- 
cation error. For q — 2, there is a marginal detection of 
non-zero spin at late times t > 3000M. 



We interpret the monotonically increasing spin Xi as 
evidence of tidal spin-up of non-rotating black holes. To 
investigate this process in more detail, we consider the 
spin xi as a function of the orbital frequency. Alvi [95] 
derived tidal spin-up as a function of binary coordinate 
separation b/M . Converting his formula into a function 
of the orbital frequency (which heuristically should be 



Xi ~ Xi, 



T]M_ 



4M 



r(i + 3xi,oo) 



(30) 



Xi, 



\ 1 ( (Mfi) 4/3 



2ri )OC 
7M 



(MO) 



7/3 



X1.00 is the spin magnitude of black hole 1 at infinite 
separation, and n )00 = M\{\ + Jl — x\,oo) ^ s tne cor " 
responding horizon radius. Dropping terms quadratic in 
Xi,oo because of their small size, this equation simplifies 
to' 



\i = \i.xl1-^(mn)V8 



+ v S imn)7/3 - (3i) 



Furthermore, the expression in parentheses in the first 
term on the right hand side is so close to unity that the 
deviation from unity is irrelevant given the small value 
of Xi,<x>- Approximating this parenthesis by unity, we 
finally find 



Xi = Xi,oo + A (Mil) 



7/3 



with the coefficient 



h 



7M 2 



7(1 + q) 4 



(32) 



(33) 



Therefore, we see that the spin xi(Mil) should follow a 
power law in frequency Mil. 

The magnitude of the change in the spin is determined 
by the coefficient fi(q), which is plotted in Fig. 12 The 
red circles denote the values of this coefficient for the 
large black hole in our simulations: The mass ratios con- 
sidered here all result in almost maximal tidal coupling, 
for maximal spin-up of the large black hole. In contrast, 
the black crosses denote the spin coupling coefficient for 
the small black hole. The spin coupling coefficient for the 
small black hole is smaller by a factor between 4 (q = 2) 
and 36 (q = 6), indicating that the smaller black hole 
will be much less susceptible to tidal spin-up. Therefore, 
from the perturbative analysis of tidal coupling, we ex- 
pect that the larger black hole in all our simulations will 
be spun up by approximately similar amounts, and that 
the small black hole will be spun up significantly less. 
This expectation is already borne out in Fig. |11[ where 
we were able to resolve the spin- up of BH 1, but not the 
(smaller) spin- up of BH 2. 

Fitting the numerical data Xi(Affi) to the functional 
form of Eq. (32 1 with the one free fitting parameter Xi,oo 
results in a moderately good fit. The fit can be im- 
proved if the coefficient /1 is also fitted for, and can be 
improved further by also allowing the exponent to vary, 
i.e. a power-law fit with an offset. The results of these 
fits (which we refer to as Fit 3, Fit 2, and Fit 1, respec- 
tively), 



are shown in Table III 



Figure 13 plots the fits 
and their residuals for mass ratios q — 2 and q = 6. All 
fits were performed over the numerical data up to orbital 
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Xj (more massive black hole) % 2 (l ess massive black hole) 




2000 4000 2000 4000 2000 4000 2000 4000 

t/M t/M 



FIG. 11. Convergence test of the dimensionless black hole spins % = S/M 2 . The left panel shows data for the more massive 
black hole, the right panel for the less massive black hole. For each mass ratio, three resolutions are shown, labeled (N=3,4,5). 
The spin of the more massive black hole, xi is convergently resolved and is monotonically growing during the simulation. The 
spin of the smaller black hole X2, is consistent with \2 = within numerical errors. 



0.01 



0.001 r 




O q=2,3 t 4,6 

x q=l/2, 1/3, 1/4, 1/6 



=M 1 /M 2 



FIG. 12. The coupling coefficient / that determines the 
magnitude of the change of the spin xi during the inspiral 
as a function of the mass ratio q = M1/M2. The red circles 
denote the coefficients for the large black hole for the mass 
ratios simulated here. The crosses denote the coefficient for 
the small black hole, which can be obtained from the same 
plot at the inverse mass ratio. 



frequency MCI — 0.05fj^] As can be seen from the in- 
sets of Fig. [13J the more general Fit 1 is superior to a fit 
with fixed exponent 7/3 (Fit 2), which in turn is supe- 
rior to the one-parameter Fit 3 of Eq. (32 1. For q = 2, 



the residual of Fit 1 is almost two orders of magnitude 
smaller than for fits 2 and 3. Coefficient A2 in Table [TTT] 
shows that the numerical data prefers a power law with 
a slightly larger exponent of roughly 8/3 instead of the 
expected 7/3. If the exponent is fixed to 7/3, then coeffi- 
cient Bi indicates that the overall magnitude of the spin- 
evolution is larger in the numerical simulation by about 
a factor of 1.3 relative to the expected behavior Eq. (32 1. 



All fits indicate fairly consistently that the spin of the 
large black hole at infinite separation would be around 





Fit 1 


Fit 2 




Fit 3 


q 


A + Ai(MQ) A2 


B + Bi/i(Mfi) 7/3 


C0 + MMQ) 7 / 3 




W 6 A A! A 2 


10 6 B 


Bi 


10 6 C 


2 


-0.95 0.0362 2.57 


-1.26 


1.23 


-0.88 


3 


-1.10 0.0496 2.64 


-1.59 


1.29 


-0.99 


4 


-1.25 0.0474 2.62 


-1.78 


1.34 


-0.96 


6 


-0.81 0.0602 2.74 


-1.39 


1.34 


-0.74 



TABLE III. Fitting parameters for fits to the xi(MCl) data. 



Beyond this frequency, we modify the gauge in the simulation, 
which leads to artifacts in xi(^A^)- 
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MQ. 




MQ. 



FIG. 13. Dimensionless spin xi of the larger black hole as 
a function of the orbital frequency MQ. Plotted is the nu- 
merical data, and three fits to the data, fitted in the interval 
MQ < 0.055. 



FIG. 14. Mf/M as function of q. Also shown are the results 
from the fitting formula of Tichy and Marronetti [96] , the an- 
alytical prediction of Berti et al. [57] , and the fit of Buonanno 
et al. [T7J to numerical data. 

mine that Einstein's equations are solved with sufficient 
accuracy. This we have done. However, in addition to 
this numerical convergence test, the resolution of the ap- 
parent horizon finder must be varied to ascertain that 
the apparent horizon is found with adequate accuracy. 
And finally, the resolution of the eigenvalue solver that 
computes the approximate helical Killing vectors on the 
apparent horizon (cf. Appendix of |81| ) must be varied 
to check that the approximate Killing vectors are cal- 
culated accurately enough. Unfortunately, we did not 
output enough data during the numerical evolutions to 
perform the second two convergence tests. 

In addition, further work would be needed to ascertain 
that the approximate Killing vectors (and the spin com- 
puted using these, cf. [5T]) are indeed generating a spin 
compatible with the spin definitions of the perturbative 
work [95j . Because of all these cautionary comments, and 
insufficient numerical data, we postpone quantitative re- 
sults about tidal spin-up to future work. 



F. Remnant properties 



10 6 , anti-aligned with the orbital angular momentum 
(cf. coefficients A a ,B ,Co). 

These results are enticing and suggestive. However, 
we caution the reader that the observed effects are very 
small, with changes to the dimensionless spin of order 
10~ 5 . Before drawing firm conclusions, one must estab- 
lish that the numerical data is accurate enough by per- 
forming a three-fold convergence test. First, the resolu- 
tion of the numerical evolution must be varied to deter- 



Figures 14 and 15 show the mass and spin of the rem- 
nant black hole (computed using approximate Killing 
vectors on the apparent horizon 81, 84-86]) as a function 
of mass ratio q. These quantities are also listed in Ta- 
ble [I] Several fitting formulas in the literature give good 
agreement with the remnant spin and are plotted in Fig- 
ure [TBJ Analytical predictions of the final mass do not 



agree as quite as well, as seen in Figure 14 however, the 



16 



0.7 



S 06 
'5. 

GO 

1 0.5 



0.4 



©— © Our simulation 
X -X Barausse, Rezzolla 09 
D d BKL 

Tichy, Marronetti 




.01 

o 

-0.01 

■0.02 









> 



0.1 



0.15 



0.2 



0.25 



FIG. 15. Sf/Mf as function of symmetric mass ratio r\. 
Also shown are the results of fitting formulas and estimates 
from Barausse and Rezzolla [98] . from Buonanno, Kidder and 
Lehner (BKL) [99], and from Tichy and Marronetti [96]. The 
inset shows the difference Xflt — Xnr between fitting formula 
and our numerical results. 



method would require smoothing put in by hand. 

The recoil speed Wkick = \v\ of the remnant is listed in 
the last column of Table [H We estimate several sources 
of uncertainty, which are listed in Table |IV| Numerical 
truncation error is estimated by taking the difference of 
«kick computed using the highest and second-highest nu- 
merical resolutions; this is the dominant source of error 
for two of our simulations. The uncertainty in extrapo- 
lating the waveform to infinity is estimated by comparing 
Wkick computed using waves extrapolated using 3rd order 
polynomials |91j versus an identical calculation using 4th 
order polynomials. The error associated with truncating 
Yi m modes for t > 6 in the momentum flux is estimated 
by comparing with an identical calculation where we re- 
tain only I < 5. Initial data effects such as the initial 
pulse of junk radiation add a spurious recoil of about 1 
to 2 km/s, depending on the run. There is an additional 
small error that results from neglecting the recoil that 
occurs in the early inspiral between t = —00 and the 
start of our simulations; this neglected recoil can be es- 
timated to 2PN order using Eq. 22 of ref. [102) . which 
yields about 0.5 km/s for the cases shown here. Figure 16 
plots the recoil versus mass ratio for our simulations and 
for two fitting formulas in the literature. We find good 
agreement. 



IV. DISCUSSION 



formula of Buonanno et al. [17] , which is a fit to numer- 
ical relativity results, shows better agreement. 

For unequal-mass binaries, linear momentum is car- 
ried off anisotropically by gravitational waves, leading 
to a recoil of the remnant black hole. The recoil speed 
of the remnant can be computed from the gravitational- 
wave momentum flux at infinity. To do this, we start 
with the Newman-Penrose quantity ^4, extracted from 
our simulations and extrapolated to infinite radius using 
the procedure of Boyle and Mroue |91j . The momentum 
flux depends on the first time integral of ^4, and com- 
puting this time integral requires two integration con- 
stants, which we determine by the procedure outlined 
in Appendix B of Ref. 100 . This procedure involves a 
minimization over a time interval [ii,t2], where t\ and 
£2 can be chosen arbitrarily. We find that varying the 
integration-constant parameters t\ and £2 hi the range 
h e [1000M, 1400M] and t 2 € [2600M, 3000M] changes 
Vvick by only a tenth of a percent. Once we have the 
time integral of ^4, we compute the gravitational- wave 
momentum flux by the procedure of Ref. [101] . keeping 
all Yi m modes through I — 6. The time integral of the 
momentum flux gives the total radiated 3-momentum P, 
and the recoil velocity is v = —P/Mf. Note that the 
recoil velocity can alternatively be computed by a time 
derivative of the Regge-Wheeler-Zerilli strain hg m rather 
than a time integral of ^'4. We use the latter method 
because differentiation amplifies noise in the waveform 
to the extent that for the runs shown here, the former 



This paper accomplishes several tasks with regard to 
simulations of BBH systems. Section |IIB| introduces an 
efficient formalism to perform root-finding necessary to 
achieve desired initial data parameters (masses, spins, 
center-of-mass frame). Each function evaluation during 
root-finding is an entire (expensive) initial-data solve, so 
it is imperative to be able to perform this procedure with 
as few function evaluations as possible. The procedure 
introduced here, based on approximate Newton-Raphson 
iteration, performs very well. As Fig.[T]shows, one or two 
high-resolution initial data runs are sufficient. Since the 
high-resolution solutions dominate the overall CPU cost, 
root-finding can thus be accomplished with marginal ex- 
tra cost. This procedure has since then been extended to 



q 


"kick 


<^kick 


<^kick 


6v kiS 


<^'k°k 


O^kick 


2 


148 


0.7 


0.4 


1 


1 


0.4 


3 


174 


6 


0.4 


0.2 


2 


0.6 


4 


157 


1.2 


0.4 


0.3 


2 


0.6 


6 


118 


4.5 


1 


3 


1 


0.4 



TABLE IV. Recoil velocity and uncertainties in km/s. Un- 
certainties (left to right) are numerical truncation error, error 
in extrapolating waveforms to infinity, the effect of using only 
a finite number of Yt m modes to compute the momentum flux, 
error involving initial transients (e.g. junk radiation), and the 
estimated recoil accumulated from t — —00 to the start of our 
simulation. 
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FIG. 16. |«kick| as function of q. Also shown are the results 
of fitting formulas and estimates from Baker et al. |103] and 
from Gonzalez et al. 14 . 



superposed Kerr-Schild data [51] . 

We then give technical details about how to simulate 
unequal-mass binaries with multi-domain spectral meth- 
ods. In particular, we extend the dual-frame formalism 
and control systems to unequal masses, introduce eccen- 
tricity removal for unequal-mass binaries, and describe 
algorithmic modifications performed during merger and 
ringdown. 

The largest part of this paper documents a new se- 
ries of unequal-mass, non-spinning BBH simulations with 
mass ratios q — 2,3,4 and 6, lasting between 15 and 22 
orbits before merger. We show that these simulations 
have high accuracy, comparable to that of the equal- 
mass simulation presented in [42| [52]. The total mass 
is conserved during the inspiral to a few parts in 10 6 (cf. 
Fig. [5]) , a convergence test on the (not time-shifted grav- 
itational wave phase) indicates that errors in our second 
highest resolution run are a few tenths of a radian. Given 
how much more challenging a mass-ratio 6 simulation is, 
we are very encouraged that the errors are only larger 
by a factor of 4 relative to the equal-mass simulation, cf. 
Fig. [6] By moving the outer boundary, we establish fur- 
thermore, that effects due to the outer boundary arise at 
the smaller level of ~ O.Olrad in the waveform, as shown 
in Fig. [8] We also perform a convergence study on the 
subdominant (3,3) and (2,1) modes of the gravitational 
radiation. These subdominant modes become more im- 



portant with higher mass ratio (see [3TJG2] and Fig. 10), 
and we argue that this increases the need for reflection 
minimizing boundary conditions, as those applied here. 
The final waveforms, extrapolated to infinite extraction 
radius, are shown in Fig. [9] 



We then consider carefully the change in the spin of the 
larger black hole. This change is broadly consistent with 
perturbative calculations of black holes: The power law 
of the spin vs. orbital frequency is rather well matched 
(~ 2.66 vs. 7/3), and the amplitude of the change is also 
reasonably close, being off by a factor ~ 1.3. A more de- 
tailed comparison must, however, await more complete 
convergence data, to allow comprehensive quantification 
of the error in the numerical spin. But nevertheless, these 
data point to the fact that our simulations are in fact for 
a BBH where the larger black hole started at infinite 
separation with a spin of ~ 10~ 6 anti-aligned to the or- 
bital momentum. Tidal spin-up increases this spin during 
the early (not modeled) inspiral, so that the spin passes 
through zero when our simulations commence. 

Finally, we compare remnant properties and kick ve- 
locities. These are found to be in reasonable agreement 
to various fitting formulae in the literature. 

An important result of this work is the accurate cal- 
culation of long subdominant mode waveforms. These 
are needed for parameter estimation, calculating physical 
quantities such as the gravitational recoil, and for mod- 
eling analytic and phenomenological waveforms (see [21] 
and references therein). Furthermore, recent results in- 
dicate that they are important for LIGO event detec- 
tion: Brown, Kumar and Nitz (in prep 2012) have found 
that for q > 1.8, the top subdominant modes must be 
taken into account in order to achieve the usual signal 
to noise ratio loss criterion "overlap greater than 0.965" . 
Pertinent factors used in these simulations which have 
contributed to the achieved accuracy are: (i) our use of 
semi-absorbing boundary conditions combined with the 
location of the outer boundary, (ii) extrapolation to infin- 
ity, (iii) good numerical resolution because of the length 
scale problem (which becomes more severe for the sub- 
dominant modes), and (iv) pseudo-spectral methods. In 
sum, we have been able to perform the first long and 
accurate numerical simulations of unequal non-spinning 
binary black holes with mass ratios as high as 6, with ex- 
cellent convergence and modest computational cost, even 
for the subdominant modes. 
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Appendix: Non-overlapping spectral grid 

In our spectral evolution code, the use of overlapping 
grids sometimes leads to weak instabilities. We find that 
these instabilities can be cured by use of non-overlapping 
grids. There are a number of choices one has to make 
while designing such a grid. A basic assumption is that 
at some distance from the center the geometry of the 
spacetime is close to spherical symmetry. Spherical shells 
are our most efficient grid structure to represent such a 
region. In the near zone (around each singularity) we 
have an excision boundary of topology S 2 which sug- 
gests that, at least in the neighborhood of each excision 

Let 



boundary, one can use spherical shells (see Fig. 17 



Ra and Rb be the outer radius of the region around the 
excision boundaries that is described by spherical shells. 
And let the coordinate centers of the excision bound- 
aries, as set by our initial data solver, be {xa, TJa, %a) and 
[xbi Vb,zb)- Assume for the simplicity of the discussion 
that xa > xb and \xa\ < \ x b\- We center the outer 
shells at the origin of our coordinate system. The inner 
radius Rc of the outer spherical region is set to approx- 
imately three times the distance between the centers of 
the excision spheres. Next we need to fill in the space be- 
tween the outer sphere <S^[(0, 0, 0), Rc] and the two inner 
spheres S a [(xa, Ua, z a ), Ra] and S B [(x B , Vb, z B ), Rb}- 

In order to construct the actual subdomains filling up 
the space between S A , S B and Sq, we will make use of 
(6, 4>) coordinates aligned with the x axis, defined with 
respect to the centers of either S EA or S EB (these spheres 
will be defined below): 



6 A = tan 1 (z/y), 



x - x E a 



yj{x- X E a) 2 + V 2 + Z 2 



(A.l) 
(A.2) 



with similar definitions for {4>b,^b)- 

We next define a projection map used to connect var- 
ious surfaces with spheres (see Fig. 18). Let Sl be a 
surface parametrized by (9,4>). Let Sfj be a sphere, and 
let Pw be a point in the interior of the sphere but not 
on the surface, Pw £ Sl- 

For each point Qu(x L ) G Sl we construct a line con- 
necting Ql and Pw- This will intersect the sphere in 
two points. Let Qu{x\j) be the intersection point that 




FIG. 17. Schematic geometry of the spherical regions of 
the grid geometry. The outer radii of the regions around the 
excision boundaries covered by spherical grid is indicated by 
arrows. The excision boundaries themselves are marked by 
circles drawn with dashed line. 




FIG. 18. Schematic diagram of the projection map used to 
connect various surfaces with spheres. Given a a surface Sl, 
a point Pw and a sphere Sfj, the projection of a point 
Ql £ cSl is defined by the intersection of the line crossing 
Pw, Ql and the sphere Sfj, such that Ql is between Pw and 



is on the same side of Pw as Ql- Thus we have de- 
fined a rule that associates a unique point Qu € Sy to 
each point Q L € Su- We will label the point Qu by the 
same parameters (#,</>) as the associated point Qu- The 
projection map is defined as 



M(P w ,Su) 



(A.3) 



P 



4(0,0) 



where we used p as a radial parameter, with range p € 
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FIG. 19. Left: Schematic diagram indicating the 

x =const plane separating the regions around the two ex- 
cision boundaries. Right: Schematic diagram of the spheres 

53 o3 q3 
BA) SeEi&EB- 



■1,1]. We have 



M(-i,e,cf>) = xi(e,< 
M(+i,e,(f>) =x\j{d,< 



(A.4) 
(A.5) 



We associate one projection map with each of the three 
spheres: 



M c 
Ma 
Mb 



= M((x c ,0,0),S 3 c ) 
= M{x\,S A ) 
= M{x%,S%) 



(A.6) 
(A.7) 
(A.8) 



whe re xg is defined in Eq. (|A.9|. As pointed out in 
Sec. 



II B 



A/B 



are slightly offset from the x axis along 
the y direction. 

Next we divide the volume in the interior of S c , outside 
of S A and into wedges of various shapes. First we pick 
an x =const plane, Vc (see Fig. 19), that separates the 
regions around the two excision boundaries, using 



x c = r?(l - x a + rjixs-, with 



XA 



4 \x A \ + \x B \ 



(A.9) 
(A.10) 



Our preferred value for n is 0.99. 

When £ < 1/3 (corresponding to mass ratios q < 2) 
we start by constructing a sphere <S| A [(^ea, 0, 0), -Rea] 
with 

x EA = 0.9nx A (A.ll) 
Rea = \/(xea - x c ) 2 + {-qxA - x c ) 2 - (A.12) 
The sphere 5j| A intersects the plane Vc in a circle 



s me '■- $ea n Vc 



with radius 



''ME 



\r]x A - x c \ 



(A.13) 



(A.14) 



On the other side of Vc we define two concen- 
tric spheres (see Fig 19): [(£eb, 0, 0), -Reb] and 




FIG. 20. Schematic geometry of the touching grid geometry. 
In a typical simulation we surround the shown grid geometry 
by about 20 further spherical shells on the outside, which are 
not shown in this diagram. 



See [(zeb, 0, 0), -R E e] with 



XEB 


= r\x B 




(A.15) 




= ?'me x max ^0.4, 


f] x B - x c 
1]X A - x c 


) (A.16) 


Ree 


= yj (xeb - x c ) 2 + 


r 2 
'me 


(A.17) 


Reb 


= J (xeb - x c ) 2 + 


r 2 

'MB' 


(A.18) 



These choices imply that £>| B intersects Vc in a circle 
with radius tmb 



< 2 



MB 



(A.19) 



Next we define wedges/cylinders filling up the space be- 
tween the three spherical surfaces. (In our terminology 
wedges have topology I 1 x B 2 , cylinders have topology 
I 1 xS 1 xI 1 .) 

• we connect the x > xc + (3/2) (xea — xc) portion 
of <Sf A with S c using Mc and call this the CA 
wedge 

• we connect the same portion of S EA with S A using 
Ma and call this the EA wedge 

• we connect the xc < x < xc + (3/2)(ieea — xc) 
portion of S EA with S c using Mc and call this the 
CA cylinder 

• we connect the same portion of S EA with S A using 
Ma and call this the EA cylinder 
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• we connect the points x l £ Vc inside S^ AE but out- 
side S^ AB with S A using A4a and call this the ME 
cylinder 

• we connect the points x % € Vc inside S^b with S A 
using M.a and call this the MA wedge 

• we connect the same set of points with S B using 
M.B and call this the MB wedge 

• we connect the x < xq — (3/2)|a; EB — x c\ portion 
of S EB with S c using A4c- The portion inside S EE 
is the EE wedge, the portion between S EE and 5^, 
is the CB wedge. 

• we connect the same portion of S EB with S B using 
Mb and call this the EB wedge. 

• we connect the xq > £ > xc — (3/2)|xeb — xc\ 
portion of S EB with S c using Mc- The portion 
inside S EE is the -EE cylinder, the portion between 
S EE and S c is the CB cylinder. 

• we connect the same portion of <S EB with S B using 
M.B and call this the EB cylinder. 

In the cases where £ > 1/3 (corresponding to mass 
ratios q > 2) we use a slightly simpler algorithm: we 
start by constructing 6> EB [(x EB , 0, 0), i? EB ] with 

x E B = r]X B (A.20) 
R EB = y/2 X |:e E b - x c \ ■ (A.21) 



The sphere S EB intersects Vc in a circle 

Smb := Sl B n P c (A.22) 

with radius 

7-mb = -xc\. (A. 23) 

On the other side of Vc we define 

xea = Vxa (A. 24) 

Rea = ^(x E A-x c ) 2 +rl lB . (A.25) 

Once again, 5 EA [(xea, 0, 0), J?ea] intersects Vc in a cir- 
cle 

S^ B :=S EA nV C . (A.26) 

The definition of the various wedges and cylinders in 
this case is similar to what is used for £ < 1/3 with the 
exception that there are no EE or ME cylinders/wedges, 
asSl, A nVc = S EB nVc=Si lB . 

See Fig. ^ for a 3D snapshot of a grid used for a run 
with mass ratio 2. This simulation uses the more com- 
plicated domain decomposition, although it is close to 
the dividing line £ = 1/3 where we switch to the simpler 
domain decomposition. As a last remark, in the runs de- 
scribed here, we have subdivided each wedge (of topology 
I 1 x B 2 ) into five distorted cubes. 
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